Angle Dependence of Landau Level Spectrum in Twisted Bilayer Graphene 
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In the context of the low energy effective theory, the exact Landau level spectrum of quasiparticles 
in twisted bilayer graphene with small twist angle is analytically obtained by spheroidal eigenvalues. 
We analyze the dependence of the Landau levels on the twist angle to find the points, where 
the two- fold degeneracy for twist angles is lifted in the nonzero modes and below/above which 
massive/massless fermion pictures become valid. In the perpendicular magnetic field of 10 T, the 
degeneracy is removed at ^dog ~ 3° for a few low levels, specifically, O^cg — 2.56° for the first pair of 
nonzero levels and Odeg — 3.50° for the next pair. Massive quasiparticle appears at 9 < 9c — 1.17° 
in 10 T, which match perfectly with the recent experimental results. Since our analysis is applicable 
to the cases of arbitrary constant magnetic fields, we make predictions for the same experiment 
performed in arbitrary constant magnetic fields, e.g., for _B = 40T we get 9c — 2.34° and the 
sequence of angles 6'dcg ~ 5.11, 7.01, 8.42, • ■ • for the pairs of nonzero energy levels. The symmetry 
restoration mechanism behind the massive/massless transition is conjectured to be a tunneling 
(instanton) in momentum space. 
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I. INTRODUCTION 

One of the most prominent hallmarks of the existence 
of massless Dirac fermions in monolayer graphene was the 
experimental discovery of an unusual quantum Hall effect 
In other words, the observation of the Landau level 
(LL) pattern of the massless Dirac fermions, in the pres- 
ence of a perpendicular magnetic field B, En ^ ±VBn 
(n = 0, 1, 2, • • • ), confirmed the characteristic band struc- 
ture of graphene, long after its theoretical prediction Q. 

The quasiparticles in bilayer graphene (BLG) of Bernal 
stacking are massive and described by the following 
Hamiltonian. For a given valley, say K, in the low energy 
continuum limit, it is 
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^K- —\g2 



9=^ 



(1) 



where m is the effective mass of the quasiparticle in BLG 
and we have introduced complex coordinates z = {x + 
iy) I \/2 on the graphene plane. The derivative operators 
are defined as 9 = djdz and d = d/dz, respectively. 
The LL spectrum of BLG is thus different from that of 
monolayer graphene and is given by eigenvalues of the 
following Hamiltonian in a magnetic field B: 
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-huj 





i2 



(2) 



that is, by En = zthtj ^/n(Jl^^^^T) Q, where w = eB/m is 
the cyclotron frequency. The lowering/raising operators 
a/a^ satisfy the usual harmonic oscillator algebra and are 
given by suitable combinations of covariant derivatives 
representing the magnetic field B. 
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When graphene layers stack together, the interlayer 
couplings significantly change the nature of quasiparti- 
cles, like the BLG case considered above. Surprisingly, 
massless Dirac fermions survive the stacking in multi- 
layer structure grown on SiC Q- The main reason for 
this stunning effect is thought to be the decoupling of 
twisted layers @ , @- 14 1 . Persistence and properties of 
massless Dirac fermions at small twist angles are, how- 
ever, under some controversy. According to ab initio 
calculations the decoupling occurs at any values of 
twist angle and massless Dirac fermions are essentially 
those of monolayer graphene. On the other hand, the 
tight-binding analysis [5| indicates a strong role played 
by the interlayer coupling, which considerably affects the 
nature of quasiparticles. These results are based on the 
band structure calculation in the absence of magnetic 
field. In a recent experiment [l6j . the issue of angle de- 
pendence of the LL's in the presence of magnetic field is 
addressed by combining scanning tunneling microscopy 
and LL scanning tunneling spectroscopy. They measured 
some critical angles at which two-fold degeneracy due to 
the presence of two massless Dirac fermions is lifted and 
where the massless Dirac fermion picture breaks down. 
Specifically, the degeneracy can be seen at angles above 
roughly 6'dcg ~ 3° for a few low levels in the magnetic 
field of 10 T, and the critical value of twist angle be- 
low/above which massive/massless LL spectrum is shown 
is dc — 1.16° for 10 T. An effective Hamiltonian to obtain 
the corresponding LL's is suggested recently but only the 
zero-modes are analytically investigated [17] . In order to 
study the angle dependence we should have controls over 
the nonzero-modes. 

In this paper, we exactly solve the LL spectrum of the 
Hamiltonian proposed in 17] to give the concrete values 
of angles, which show very precise agreements with the 
measured values for the magnetic field 10 T [l6|. Since 
the LL spectrum is obtained in analytic form, we can 
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predict 9c and 0dcg for every nonzero pair of the LL's in 
the presence of arbitrary constant magnetic field (under a 
plausible assumption on the physical continuity of the LL 
spectrum as functions of twist angle about which we dis- 
cuss in Appendix IXj) . Furthermore, our analytic result 
for LL's (fT8|) smoothly interpolates between the spec- 
tra known before for massive/massless quasiparticles in 
BLG, and allows one to get systematic power series cor- 
rections for both of two sides of the spectrum. 

A natural question one can ask here is: "What is the 
symmetry restoration mechanism behind the transition 
between the massive/massless spectra?" Based on our 
exact results, we anticipate that the non-perturbative 
symmetry restoration mechanism is a tunneling (instan- 
ton) in the momentum (reciprocal) space [Tsj . 

The remaining parts of this paper are organized as fol- 
lows. In Sec. nil we briefly review the construction of 
the low-energy effective Hamiltonian for charged quasi- 
particles in twisted BLG [l3| and solve the associated LL 
problem analytically by invoking a differential-equation 
representation of the eigenproblem. In Sec. IIIIl we re- 
veal the distinction between two asymptotic regions in 
LL spectrum, obtained in Sec. [TTl considered as func- 
tions of twist angle. Finally, in Sec. IIVI we wrap this 
paper up with a short summary and some discussions. 
An appendix is devoted to discuss the change in the LL 
spectrum driven by the Fermi speed renormalization as 
twist angle increases and to advocate the validity of the 
analysis made in Sec. IIIIl 



II. EXACT LANDAU LEVELS IN TWISTED 
BILAYER GRAPHENE 

Twisted bilayer graphene is a structure specified by a 
rotational mismatch given by an angle 9 with respect to 
the perfect Bernal (AB) stacked bilayer graphene. This 
twisted pattern is not difficult to find but can be seen 
on the surface of graphite, for an example. In low en- 
ergy limit, it shows a drastically different electronic band 
structure from that of the Bernal-stacked BLG. Its low 
energy quasiparticles are two massless Dirac ferniions 
rather than one massive fermion, per each valley (K/— K) 
0. For a small 9, the apices of the associated Dirac cones 
are separated by |AK| |K — Ke| ~ 47r0/3-\/3acc in 
reciprocal space, where Occ — 1.42 A is the distance be- 
tween two adjacent carbon atoms in the hexagonal lat- 
tice. A commensurate rotation with a periodic Moire 
pattern occurs at the angles 9i: 



The reciprocal lattice is spanned by 



cos 9i 



ii^ + 3i + 



ii^ + 2,i + l 



0,1, 2,--- , 



(3) 



and the superlattice structure is specified by basis vectors 
ti — iai-t-(i-|-l)a2 andt2 = — (i-f l)ai-|-(2i-|-l)a2, where 
ai, a2 are the Bravais lattice basis vectors in the single 
layer hexagonal lattice [s^. The lattice constant of the 
superlattice is given by |ti| = |t2| = Occ ^9^2 -f 9i + 3 . 



Gi 



47r 



-[(3i + l)ai +a2] 



9^2 -(- 9i -I- 3 ' 

G2- Q^2^^9, + 3 [-(3» + 2)ai + (3» + l)a2]. (4) 

The first Brillouin zone for twisted BLG is depicted in 
Fig. [iKa) and the corresponding low energy band struc- 
ture in the K- valley is shown in Fig. [Ub). Electronic 
properties of twisted BLG and related systems are cur- 
rently under intensive debates l6i-ll4ll. 




FIG. 1: (Color online) (a) First Brillouin zone for twisted 
BLG. The first Brillouin zone of the upper layer (dashed 
hexagon) is rotated by an angle 6 with respect to that of the 
lower layer (full hexagon), (b) Low energy band structure 
near the K-valley of twisted BLG. The dispersion relation for 
quasiparticles is given by E{k, k) = zfc J^^^^ \ k — ^^\\k + ^^\, 
where k = -^{k^+iky), AK = -^{AK^ + iAKy), and m{e) 
is 0-dependent effective mass. If we set k = k — then we 
get the massless behavior E ~ ±\/2/i{)p|K| near k — 0, where 
vf = is the renormalized Fermi speed. 

\/2m(y) 

In a recent paper [l^l , an effective Hamiltonian for low 
energy dynamics of twisted bilayer graphene is proposed. 
Neglecting commensuration effects between two layers, 
the Hamiltonian describing twisted BLG around the pair 
(K,Ke) reads 



/^tw(k) = 



Hoik 



(5) 



where Hd is the Dirac Hamiltonian for monolayer 
graphene and H± is the hopping matrix between the two 
layers. According to the analysis of the Moire pattern in 
twisted BLG 0, for small twist angle 9, there are three 
different types of H±: 



1 1 
1 1 



e**'3 e 



1 



(6) 



where ij_ is a coupling parameter which generally de- 
pends on 9. The 6'-dependence of tj, is, however, very 
mild and the Slater-Koster calculation performed in 
Ref. indicates that f_L — 0.4 71 is nearly a constant, 
where 71 ~ 0.3 eV is the nearest interlayer coupling in 
the untwisted Bernal-stacked BLG. The first one in ^ 
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corresponds to the Fourier component with the crystal 
momentum G = and the second and the third one 
with G = — Gi, — (Gi + G2). The summation in the in- 
terlayer coupling term in ([5|) runs over these three types 
of the coupling matrices and the other components are 
suppressed. Due to this interlayer coupling, there is a 
slight difference of 1 meV order between the energies as- 
sociated with each Dirac points but this second-order ef- 
fect in perturbation is negligible. Assuming a simplified 
interlayer coupling under the condition t± ^ fiwF|AK| 
{vp is the Fermi speed in the monolayer graphene) [2^ , 
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(7) 



one can obtain an effective two-band Hamiltonian which 
resembles (P) 



rroff 







m{9) ~ (A^/2) 



where m{9) — 15i±/Avp is 0-dependent effective mass 
due to the 0-dependent t±. In ^ the multiplication 
by 3 mimics the summation over three possible crystal 
momentum G, and another factor 5/2 is introduced to 
match the spectrum at 6* = (Bernal-stacked BLG). For 
6 = 0, the period of the superlattice is infinite and the 
summation over G is overcounting since Gi = G2 = 0. 
Thus the multiplication by 3 in ([T]) should be disregarded 
in this case and the interlayer coupling becomes that of 
the Bernal-stacked BLG: 



Hde = o) = -^ 



1 




-71 




(9) 



Then the effective Hamiltonian goes to ([1]). Since 
|AK| is proportional to twist angle 6, the two-band ap- 
proximation made here is restricted to be applicable for 
very small angles only. Nevertheless, we believe that the 
LL spectrum given by (1181) below is smoothly connected 
to the spectrum for larger angles, and the analysis made 
in Sec. |lll] is trustable. The multiplicative factor in ([7]) 
also plays a crucial role in this context. We relegate the 
discussion of this issue to the appendix of this paper. 

In the presence of a perpendicular magnetic field, B, 
the Hamiltonian is written in the form 



Hp 



-huj{9) 
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t2 



(10) 



where uj{0) = eB/m{6). The lowering/raising operators 
a/a^ are given in terms of the covariant derivatives D = 
a+ f z) / D = d+ fA{z, z), and A= ^{A^- lAy) 
is a complex vector potential. With the gauge choice 
A = —^Bz, we specifically have 



which satisfy [a, a] = [a^,a^] = and [a, a'''] 
complex parameter proportional to /S.K as 



AeB 



■AK 



ieB 



{AK^ + iAKy 



1. /? is a 



(12) 



We can make /3 real- valued simply by rotating the coor- 
dinates, 



13 



ieB 



lAKI 



(13) 



After the rotation, for ip = {ipi 1^2)'^ , the eigenvalue prob- 
lem for the Hamiltonian ([TU| . Hb^P = Eip, reduces to the 
following one-component problem: 



(a^^-/3')(a^-/3')Vi=AV'i 



(14) 



where A = [E/huj(0)]'^. The remaining component can 
be expressed in terms of -01, A, and the lowering operator 
a (except the case A = 0). Using the anti-holomorphic 
representation [l9j . 



a I— ^ 



dx ' 



(15) 



the eigenvalue problem ([T4| is expressed as a second order 
ordinary differential equation, 



{x^ - (3^){u" - - Xi 







(16) 



Here, the variable x is not the coordinate in the graphene 
plane and therefore the function u{x) representing eigen- 
states are not the wave function in coordinate space. By 
rescaling x 1-^ (3x and setting u{x) = V x^ ~ 1 v(x), we 
get the spheroidal equation oi p = /3'^ and q = 1 



[{x^ - l)v']' + 



-/(x2 - 1) - A 



1 



v = (17) 



which is a particular case & = s = of the confluent 
Heun's equation: [{x^ — l)w']' + [—p'^{x'^ — 1) + 2pbx — X — 



SL+l^±i!l^]v = [2C 



The eigenvalues of the confluent Heun's equations are 
A^°] n-i{Pj b), and then the corresponding LL's are given 
by the spheroidal eigenvalues 



(fi^w)'Ai':L„l(/3^o), 



(18) 



where the superscript (a) stands for 'angular'. The LL 
spectrum for a fixed value of /3 (twist angle) is depicted 
in Fig. [21 This analytic result reproduces the numerical 
calculation performed in Ref. |17], except the scale of in- 
dependence due to the multiplicative factor in ([7]). 



III. ANALYSIS OF SPECTRUM IN TWO 
ASYMPTOTIC REGIONS 
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Let us consider the region of the parameter p in which 
the band structure is well-described by two massless 
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FIG. 2: (Color online) Landau levels of twisted bilayer 
graphene, given by the spheroidal eigenvalues. The character- 
istic energy scale of the van-Hove singularity has been chosen 

0.1 eV. This scale corresponds to 



as = r|AK|78m = 

the twist angle of ~ 3.27' 



Dirac fermions whose Fermi speed is renormalized to be 
smaller than that of the monolayer graphene. In this 
region, the LL's must show a two-fold degeneracy that 
reflects the existence of two copies of fermions. Fig. [2] 
reveals its tail, and Fig. [3] clearly indicates that degen- 
eracy. The lifting of this degeneracy signals breakdown 
of the description of electronic bands by two copies of 
massless fermions and it is caused by the contribution 
from a saddle point in the band located between the two 
Dirac cones. Twist-induced van-Hove singularity eventu- 
ally dominates the spectrum. 
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FIG. 3: (Color online) Landau level energies divided by Q- 
dependent energy scale fiijj{6), for a fixed perpendicular mag- 
netic field B = 10 T, plotted against twist angle. If the 9- 
dependence of ijj{6) is mild enough, the qualitative behavior 
of the LL's is the same as E/fkj{6). 



Small /3 expansion [2l| for the eigenvalues gives, 
2n{n- 1) 



El 



n{n — 1) — 



{hjj(e)Y " ^' (2n + l)(2n-3) 

2n{n - l)[n{n - l)(4n(n - 1) - 39) + 63] 



/3« 



(2?i + 3)(2n l)3(2n - 3)3(2n - 5) 
+ 0(/3''). (19) 

Untwisted (/3 — 0) LL spectrum is given by the leading 
term in (jl9p . The common factor n{n — 1) appears in 
all the higher order terms in and the modes labeled 
by n = and n — 1 remain to be zero-energy modes 
even if the corrections in higher powers of /? are consid- 
ered. Therefore they are protected. Notice that we are 
neglecting any other potential lifts of degeneracy due to 
Zeeman effect, interactions, etc. This protection is due 
to the topological structure of the band in twisted BLG 
[l7| . In twisted BLG, the two Dirac cones are not re- 
lated by time-reversal symmetry and they are described 
by the same Berry phase, so that they are of different 
topological structure from the Dirac cones in monolayer 
graphene. 

For sufBciently large /3, the asymptotic expansion of 
the eigenvalues I21| gives 
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^ + -- - , (20) 



where the square bracket denotes the integer part of n/2. 
The presence of this integer- valued function implies that 
there are asymptotic degeneracies between the each even 
level and the next odd one. The difference between the 
2f-th and [21+ l)-th LL's decreases exponentially and is 
vanishing at infinite /3: 



p2 

^21+1 



^2e 



4i+2 



{£ -iy.il 



.-2/3^ 



(21) 



The exponential factor is independent of i but the power 
of the monomial function is proportional to i so that, 
in physics experiments, the degeneracy between the ad- 
jacent two levels of lower £ (or n) is lifted for smaller /3 
(smaller 9), which is consistent with Fig.[3l From now on, 
we will set the external magnetic field i? = lOT, perpen- 
dicular to the plane of BLG. Since the value of /3 is large 
enough up to a very small value of 6* (/3 ^ 1 for 9 ^ 1.17°), 
we can use the asymptotic formula pT|) in order to esti- 
mate the point on which the two-fold degeneracy is lifted, 
under the assumption that the 0-dependence of oj{9) (or 
that of t±) is mild enough. The values /3deg for £ = 1 
where the n ~ 2 and n = 3 levels become non-degenerate 
and for £ = 2 (n = 4, 5) are, from (PTjl . 



o4£+3 fl4£+2 
^ Pd 



Meg 



{£~l)\£\ 



and thus 



1=1) 



2.18 and ^, 



1 



(fc2) 



2.99 



(22) 



(23) 



5 



respectively. They correspond to the twisted angles 
0(jeg — 2.56° and ~ 3.50° (see the Fig. [3]), in agreement 
with the measured value, about 3°, from Ref. 16]. 

As the twist angle 9 decreases, f3 crosses the transi- 
tion point /3c = 1 between the region where the small 
P expansion p9p can be trusted and the region where 
the large [3 expansion pO|) is trustworthy. We will call 
the region 13 > j3c massless region and (3 < /3c massive 
region, respectively, because of an obvious reason from 
the behaviors of the LL's in each domain, (IT^ and (I^Ul) . 
The critical point /3c = 1 can roughly be considered as 
the point at which the van-Hove singularity eventually 
starts dominating the spectrum and the description of 
the band by two massless Dirac fermions breaks down. 
This critical value /3c = 1 corresponds to the twist angle 
9c — 1.17° as mentioned already. The measured value 
^(mcasuiod) ^ 2. 16° jl^ is vcry close to our theoretical 
value, even though there is no exact criterion of fixing 
the value due to the continuity of the spectrum. 

Currently obtainable maximum value of a static mag- 
netic field is about 40 T. The critical value Pc — 1 and 
the value of /3dcg for each £ in (|22p are independent of 
B and thus the only effect on the values of various an- 
gles for the 40 T magn etic field is a multiplication by 
the factor 2 = y/W/TO from ^ to the angles for 10 T. 
Therefore 9c — 2.34° and the angles above which two-fold 
degeneracy can be seen are as tabulated below, in a high 
magnetic field i? = 40T. 
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5.11 


7.01 


8.42 





This simple dependence of the specific angles ^deg and 9c 
on the magnetic field B, that is. 



7dcg 



■ fi'dcg,c(So) 



(24) 



where Bq means a reference value of magnetic field, say, 
Bq = 10 T allows us to draw the Fig. 2] to read the angles 
off for various values of B. 



IV. SUMMARY AND DISCUSSIONS 

To summarize, we solved the eigenvalue problem for 
LL's in twisted BLG with small twist angles and ana- 
lyzed the angle dependence of the spectrum in a fixed 
magnetic field. The corresponding eigenvalues (jl8p are 
given by the so-called spheroidal eigenvalues modulated 
by an unknown function hLj{9) of twist angle 9, but the 
modulation does not change the spectrum seriously. As 
the results, we got the angle 9c below/above which the 
spectrum behaves as LL's for massive/massless fermions 
and also the angles 0dcg below which the two-fold degen- 
eracies due to the Dirac-point splitting are lifted. The 
specific values for {^c^'dog} for i? < 50T can be read off 
from the Fig. 21 They can be measured in principle by a 



5(T) 



10 20 30 40 50 



FIG. 4: (Color online) The dependence on B of the various 
specific angles. 



judicious experiment which combines the scanning tun- 
neling microscopy and the LL scanning tunneling spec- 
troscopy. 

The differences between each adjacent energy levels 
(1211) show non-perturbative behavior in the expansion 
parameter, 1//?^. Indeed, the exponential factor in ((2T|) 
gives us a hint about the symmetry restoration mecha- 
nism behind the transition described above. It is a typical 
signal of tunneling mechanism (instanton effect), in this 
case not in coordinate space but in momentum space [l^ . 
The situation is analogous to the well-known double- well 
potential problem in quantum mechanics. The two Dirac 
points correspond to classically degenerate ground states, 
and the van- Hove energy scale plays the role of the height 
of potential barrier between them. By the tunneling, the 
Dirac electrons are delocalized in momentum space, thus 
become localized in coordinate space. 

The deformation of the band structure of BLG trig- 
gered by twisting is the Dirac-point splitting, as depicted 
in Fig. [Ifb). The Dirac-point splitting also happens in 
untwisted Bernal-stacked BLG by the effect of external 
magnetic field which is parallel to the graphene plane 
(2^ [23} . Therefore, in a tilted magnetic field with re- 
spect to the graphene plane, the LL spectrum of un- 
twisted BLG is expected to be the same as ([T5|) . if the 
parameter /3 c>c B\\/^/B±_ represents the effect of the par- 
allel component to the spectrum [23]. In principle, the 
magnetic field component parallel to the plane, which is 
easy to get by tilting the graphene sample in an external 
magnetic field, can split the Dirac point and make the 
LL spectrum doubly-degenerate. It is intriguing to in- 
vestigate the combined effect of twisting and inclination 
in a magnetic field, since in reality graphene samples al- 
ways deviate from the perfect plane by various physical 
reasons. 
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Appendix A: Renormalized Fermi speed and angle 
dependence 

The authors of Ref. fr?] managed to get the two-band 
effective Hamiltonian (|8]) by assuming a simplified form 
of the interlayer coupling ([T]) under the condition that 
5 = ?iUF|AK|/tj^ ^ 1, as we reviewed in Sec. [TTl The 
reduction from the full Hamiltonian ([5]) to the effective 
one (HI is reminiscent of that of the Bernal-sracked BLG 
case. In fact, in addition to the simplified form assumed 
by the authors of Ref. [Tt} . we need a numerical multi- 
plication factor in ([7]) in order to connect smoothly the 
spectrum and the Fermi speed to those for larger angles 
1) obtained in Ref. @ as we shall see below. 

The condition of small ^, ^ <C 1, used to get the 
two-band effective Hamiltonian ([8]), can be understood 
and refined as follows. The minimum energy of the up- 
per band of the full Hamiltonian, approximately ii'min — 
15t_L/2, should be much larger than the van-Hove en- 
ergy of the lower band E^u = h'^\AK\'^ /8m{e) = 
AKp/30 so that they cannot feel the existence 
of each other (Fig. [5] (a)). This refined condition that 
i± ^ j^SijfIAKI <C 15) yields the restriction on the 
range of angles 



(a) 



(b) 



~ 45v^ac 



9.80° 



(Al) 



but this does not tell us sharply how small the twist angle 
9 should be. For example, a blind application of our 
formula to 6* = 8° gives a divergent result for a small 
but finite B, say, B = 0.1 T. Then, how can we trust 
the analysis presented in Sec. IHIP Some of the predicted 
values of ^c.deg might be in the region where the spectrum 
is not quite close to that given by p8)) . 

Let us recall that for sufficiently large angles, the spec- 
trum given by ([T5)) shows the characteristic LL behavior 
of massless particles. The combination of twist angle 
6^5° and magnetic field i? = lOT for instance corre- 
sponds to /3 ~ 4.26, in the "massless region" . Moreover, 
6^5° and /3c = 1 require B ^ 182 T, a very large field 
strength. This simple observation means that the LL 
spectrum of twisted BLG ^TE\\ with 9^5° cannot be 
pushed to "massive region" , unless a very strong mag- 
netic field is applied. Let us explain this point in more 
detail. As the twist angle increases the van-Hove energy, 
playing the role of a barrier between the two Dirac points, 
also increases as can be seen in Fig. [5](b). Therefore the 
"tunneling" between the two Dirac points is suppressed 
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FIG. 5: (Color online) Schematic plots for the (upper half) 
band structure of the full Hamiltonian ((5|) along the AK- 
direction. (a) The case in which the minimum energy of the 
upper band is larger than the van-Hove energy. The two-band 
approximation can be applied without any difficulty, (b) The 
opposite case. The condition to derive the effective two-band 
Hamiltonian is not satisfied any more. 



and the Z2 symmetry is broken. Unless a very strong 
magnetic field is applied, the broken symmetry cannot 
be restored and the spectrum is described by two mass- 
less fermions degenerate in energy, for 0^5°. This fact 
is clearly shown in Fig. [HI plotted for 6 = 5°. 




B{T) 



FIG. 6: (Color online) The LL spectrum HH)) for S = 5° . Each 
adjacent levels are degenerate and in the "massless" region if 
the applied magnetic field is not very high. 

Indeed, the masslessness of quasiparticles near one of 
the two Dirac points was shown in Ref. Q by applying 
perturbation theory under the opposite condition, that 
^ ^ 1. Their result indicates the Fermi speed renormal- 
ization 



(A2) 



where wp is the renormalized Fermi speed. This result 
also seems strange since it tells us that the renormalized 
Fermi speed vanishes at ^ = 3, corresponding to ~ 2°, 
and becomes negative below that angle though it should 
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be a positive quantity by definition. Therefore tfie spec- 
trum obtained in Ref. |5j] seems valid at most in tfie region 
^ ^ 3. At any rate, the LL spectrum according to the 
Fermi speed renormahzation shown in Ref. Q is given by 



Ei = ±vf 1-9 



'h^vl\AK\^ 



V2heBe 



(^ = 0,1,2,...) 



(A3) 



that is, the LL for massless fermions modulated by the 
renormalized Fermi speed. 

Note, however, that the validity regions for the spectra 
p8)) and (|A3[) can have some overlap for 3 ^ ^ <C 15. 
Actually, as we remarked above, both the spectra at an 
intermediate angle such as 6* = 5° are of massless charac- 
ter. Except the modulation function hui{9) depending on 
the renormalized Fermi speed, they are the same unless 
the applied magnetic field is very very large. The renor- 
malized Fermi speed according to ([S]) and (PO)) (see also 
the caption of Fig. [1]) is, to the first order in ^/15, linear 
in 



vf 



2fiuF|AK| 

15t_L 



15 



(A4) 



All these circumstances make it plausible that the spec- 
tra (fTS)) and (jA3p are smoothly connected in the inter- 
mediate range of ^ (or 9) and, here comes the punchline, 
the curves representing the renormalized Fermi speed for 
these two cases are almost tangent to each other, at the 
angle 6 ~ 3.37° (^ ~ 5.20). The linear function /(^) 
which is exactly tangent to the curve vf/vf — 1 — 9/^^ is 
/(^) = 2^/9^3. See the inset in Fig. [71 If the renormal- 
izaion of the Fermi speed is given by the curve depicted 
here, the LL's smoothly connecting (IT51) and (jA3l) should 
behave qualitatively as Fig. [T] The Fermi speed renor- 
maliztion affects the spectrum for angles 9 > 3.37° to 
change the shape of its tail. 

Recall that the values of ^^dog predicted in Sec. IIIII 
are independent of the (renormalized) Fermi speed since 
their predictions are based only on the spheroidal eigen- 
values themselves, i.e., En/fiuj. Therefore, if the discus- 
sion made in this appendix based on physical continuity 
of the spectrum is correct, the predicted values can re- 
main trustable. In other words, if one can find an exact 
formula for renormalized Fermi speed as a function of 9 
which interpolates the two asymptotic forms (jA2p and 
(|A4[) (equivalently the modulation function haj{9)), the 
exact LL spectrum is given by the spheroidal eigenval- 
ues modulated by it, ([T5)) . at any values of 9. Thus it is 



extremely interesting to find such a exact interpolating 
function for the renormalized Fermi speed. 

A final remark is in order. Since the validity of the 
asymptotic form (jA2p is uncertain in the intermediate 
region, the slope in the other asymptotic form (IA4p is also 
uncertain accordingly. Still the inclusion of the factor 
3 in ([7]) seems crucial, because otherwise the slope in 
(|A4[) must be modified into a too large number to have 
a chance of smooth interpolation. The other numerical 
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FIG. 7: (Color online) The qualitative behavior of the LL's 
(B = 10 T) as functions of 8, modulated by the renormalized 
Fermi speed. Inset: expected behavior of the renormalized 
Fermi speed, interpolating (|A2[) and (|A4[) (the solid lines). 



factor 5/2 in ([7]) is the actual source of uncertainty, and 
it can be replaced by a number in some range — roughly 
from 2 to 2.5. If we adopted the simplified interlayer 
coupling 



-3 X 2< J 



1 




(A5) 



instead of ([7]), the slope of (IA4[) becomes 1/6 which is a 
reasonable number to make the interpolation. Actually, 

the interlayer coupling matrix H± = —2i± ^ ^ ap- 
proximates each of the coupling terms in ([5]) much more 
closely, as one can see by performing the diagonalization 
of ([5]) after turning off the block diagonal Dirac Hamilto- 
nians. The coupling constant i± should be close to 7i/2 
in this case. 
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